/**
 * @file   Splines.h
 * @author ubuntu <dadi@ubuntu>
 * @date   Tue Jul 20 22:16:51 2021
 * 
 * @brief  Splines.
 * 
 * 
 */

#ifndef __Dadi__SPLINE__
#define __Dadi__SPLINE__

#include<iostream>
#include<lapacke.h>
#include<cmath>
#include<vector>
#include<string>

#define TEMPLATE template <typename T>

TEMPLATE class Spline;

TEMPLATE
class Spline
{
public:
    Spline();
    virtual T Spline_generator(T x) = 0;
};

TEMPLATE
Spline<T>::Spline()
{
};

TEMPLATE
class Complete_cubic_spline : public Spline<T>
{
private:
    int N;			/**< Number of knots. */
    std::vector<T> knot;	/**< The knots of the spline. */
    std::vector<T> m;		/**< The derivatives at N interpolation knots. */
    std::vector<T> ic;		/**< Interpolation conditions, the values at N interpolation knots. */
    std::vector<T> lambda;
    std::vector<T> mu;
    
public:
    Complete_cubic_spline();
    void set_N(int _N);
    int get_N();
    void set_knots(std::vector<T> k);
    std::vector<T> get_knots();
    void set_derivatives();
    std::vector<T> get_derivatives();
    void set_boundary_condition(T dfa, T dfb);
    std::vector<T> get_boundary_condition();
    void set_interpolation_condition(std::vector<T> IC);
    std::vector<T> get_interpolation_condition();
    void set_lambda();
    std::vector<T> get_lambda();
    void set_mu();
    std::vector<T> get_mu();
    T Spline_generator(T x) override;
};

TEMPLATE
Complete_cubic_spline<T>::Complete_cubic_spline()
{
    N = 0;
    knot = {0.0};
    m = {0.0};
    ic = {0.0};
    lambda = {0.0};
    mu = {0.0};
};

TEMPLATE
void Complete_cubic_spline<T>::set_N(int _N)
{
    N = _N;
};

TEMPLATE
int Complete_cubic_spline<T>::get_N()
{
    return N;
};

TEMPLATE
void Complete_cubic_spline<T>::set_knots(std::vector<T> k)
{
    knot = k;
};

TEMPLATE
std::vector<T> Complete_cubic_spline<T>::get_knots()
{
    return knot;
};

TEMPLATE
void Complete_cubic_spline<T>::set_derivatives()
{
    m.resize(N);
    for(int i = 0; i < N; i++)
	m[i] = 0.0;
};

TEMPLATE
std::vector<T> Complete_cubic_spline<T>::get_derivatives()
{
    return m;
};

TEMPLATE
void Complete_cubic_spline<T>::set_boundary_condition(T dfa, T dfb)
{
    m[0] = dfa;
    m[N-1] = dfb;
};

TEMPLATE
std::vector<T> Complete_cubic_spline<T>::get_boundary_condition()
{
    std::vector<T> boundary_condition(2);
    boundary_condition[0] = m[0];
    boundary_condition[1] = m[N-1];
    return boundary_condition;
};

TEMPLATE
void Complete_cubic_spline<T>::set_interpolation_condition(std::vector<T> IC)
{
    ic = IC;
};

TEMPLATE
std::vector<T> Complete_cubic_spline<T>::get_interpolation_condition()
{
    return ic;
};

TEMPLATE
void Complete_cubic_spline<T>::set_lambda()
{
    lambda.resize(N-2);
    for(int i = 0; i < N-2; i++)
	lambda[i] = (knot[i+2] - knot[i+1])/(knot[i+2] - knot[i]);
};

TEMPLATE
std::vector<T> Complete_cubic_spline<T>::get_lambda()
{
    return lambda;
};

TEMPLATE
void Complete_cubic_spline<T>::set_mu()
{
    mu.resize(N-2);
    for(int i = 0; i < N-2; i++)
	mu[i] = (knot[i+1] - knot[i])/(knot[i+2] - knot[i]);
};

TEMPLATE
std::vector<T> Complete_cubic_spline<T>::get_mu()
{
    return mu;
};

TEMPLATE
T Complete_cubic_spline<T>::Spline_generator(T x)
{
    T K[N-1];
    for(int i = 0; i < N-1; i++)
	K[i] = (ic[i+1] - ic[i])/(knot[i+1] - knot[i]);
    T b[N-2];
    for(int i = 0; i < N-2; i++)
    {
	if(i == 0)
	    b[i] = 3*mu[i]*K[i+1] + 3*lambda[i]*K[i] - lambda[i]*m[0];
	else if(i == N-3)
	    b[i] = 3*mu[i]*K[i+1] + 3*lambda[i]*K[i] - mu[i]*m[N-1];
	else
	    b[i] = 3*mu[i]*K[i+1] + 3*lambda[i]*K[i];
    }
    int M = (N-2)*(N-2);
    T A[M];
    for(int i = 0; i < M; i++)
	A[i] = 0;
    for(int i = 0; i <= N-3; i++)
    {
	if(i == 0)
	{
	    A[0] = 2;
	    A[1] = lambda[i+1];
	}
	else if(i == N-3)
	{
	    A[M-2] = mu[i-1];
	    A[M-1] = 2;
	}
	else
	{
	    A[i*(N-2)+i-1] = mu[i-1];
	    A[i*(N-2)+i] = 2;
	    A[i*(N-2)+i+1] = lambda[i+1];
	}
    }
    int ipiv[N-2];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, N-2, 1, A, N-2, ipiv, b, N-2);
    for(int i = 1; i < N-1; i++)
	m[i] = b[i-1];
    int l = 0;
    for(int i = 0; i < N-1; i++)
    {
	if(x >= knot[i] && x <= knot[i+1])
	{
	    l = i;
	    break;
	}
	else;
    }
    T p;
    p = ic[l] + (x-knot[l])*m[l] + pow(x-knot[l],2)*(K[l]-m[l])/(knot[l+1]-knot[l])
	+ pow(x-knot[l],2)*(x-knot[l+1])*(m[l]+m[l+1]-2*K[l])/pow(knot[l+1]-knot[l],2);
    return p;
};


#undef TEMPLATE

#else
//Do nothing.
#endif

